'''
Created on Nov 30, 2010

@author: oabalbin
'''

import subprocess
import snps.gatk_realignment as gra
import snps.picard_commands as mypicard
import snps.base_quality_calibration as gqc
from collections import defaultdict, deque


def lane_level_realignment(sample_lanes_dict, gatk_run_dict):
    '''
    Input: a sample lane dictionary with key:sample, values: lanes that belong to that sample
    By default this realignment method uses only known indel Sites to realign. 
    The known indel sites are provided in the indeldb_file
    gatk_run_dict is a dictionary with the parameters to run the gatk realignment walker
        indeldb_file
        
    
    '''
    # gatk parameters
    '''
    ref_genome='/exds/projects/alignment_indexes/gatk/hg19/hg19.fa'
    ref_dbsnp='/exds/projects/alignment_indexes/gatk/hg19/dbsnp132_00-All_processed5.processed.vcf'
    bam_file_name='/exds/users/oabalbin/projects/snps/exomes/aM18/test/s_3_12_sequence.hg19.aln2.rmdup.sorted.bam'
    use_mem=24
    num_cores=6
    recal_analysis_outputdir='/exds/users/oabalbin/projects/snps/exomes/aM18/test/recal_analysis/'
     
    gatk_run_dict = {'path_to_gatk':'/exds/sw/bioinfo/gatk/GenomeAnalysisTK-1.0.4705/', 
                     'resources_folder':'/exds/sw/bioinfo/gatk/GenomeAnalysisTK-1.0.4705/resources/', 
                     'rscipt_path':'/exds/sw/local/R-project/bin/Rscript'}

    '''
    
    
    # Run parameters
    use_mem, num_cores = gatk_run_dict['use_mem'], gatk_run_dict['num_cores']
    path_to_gatk, path_to_picard = gatk_run_dict['path_to_gatk'], gatk_run_dict['path_to_picard']
    ref_genome, snpdb_file, indeldb_file = gatk_run_dict['ref_genome'], gatk_run_dict['snpdb_file'], gatk_run_dict['indeldb_file'] 
    path_to_intervals = gatk_run_dict['path_to_intervals']
    recal_analysis_outputdir=gatk_run_dict['recal_analysis_outputdir']
    temp_dir = gatk_run_dict['temp_dir']
    
    new_sample_lanes_dict=defaultdict(deque)
    
    # Create realignment intervals at known indel positions
    intervals_to_realign = gra.intvTorealign_only_knownSites(ref_genome, use_mem, num_cores, 
                                                         path_to_intervals, path_to_gatk, 
                                                         indeldb_file)
    
    for sp, thlanes in sample_lanes_dict.iteritems():

        for indexed_bam_file in thlanes:
            
            realigned_bam_file = gra.realign_bam_only_knownSites(ref_genome, indexed_bam_file, 
                                                           intervals_to_realign, indeldb_file, 
                                                           use_mem, num_cores, temp_dir,path_to_gatk)
            
            dedup_bam = mypicard.markDuplicates(realigned_bam_file, use_mem, path_to_picard)
            
            recal_bam_file = gqc.main_baseQ_recalibration(ref_genome, snpdb_file, 
                                                               dedup_bam, use_mem, num_cores,
                                                               gatk_run_dict, recal_analysis_outputdir)
            
            # new dictionary with the realigned and 
            # recalibrated bam files paths for each lane that belongs to sample
            new_sample_lanes_dict[sp].append(recal_bam_file)
            
    
    for sp, nwlanes in new_sample_lanes_dict:
        # merge recalibrate bam files that belong to sample
        merged_recal_bam_file =  'merged'
        merged_dedup_bam = mypicard.markDuplicates(merged_recal_bam_file, use_mem, path_to_picard)
       
    return  merged_dedup_bam
    
    
        
'''
for each lane.bam
    realigned.bam <- realign(lane.bam) [at only known sites]
    dedup.bam <- MarkDuplicate(realigned.bam)
    recal.bam <- recal(dedup.bam)

for each sample
    recals.bam <- merged lane-level recal.bam's for sample
    dedup.bam <- MarkDuplicates(recals.bam)


for each lane.bam
    realigned.bam <- realign(lane.bam) [at only known sites]
    dedup.bam <- MarkDuplicate(realigned.bam)
    recal.bam <- recal(dedup.bam)

for each sample
    recals.bam <- merged lane-level recal.bam's for sample
    dedup.bam <- MarkDuplicates(recals.bam)
    realigned.bam <- realign(dedup.bam) [with known sites if possible]

'''
